MODULE module_wrf_getvar
  USE module_wrf_global_attr

  REAL, allocatable,dimension(:,:) :: xlat,&! latitude of model horizonal grids
                                      xlon  ! longitude of model horizonal grids
  REAL,ALLOCATABLE,DIMENSION(:,:,:):: z,&   ! height  
                                      p,&   ! pressure 
                                      ph,&  ! perturbance of pressure 
                                      phb,& ! base of pressure 
                                      pp,&  ! 
                                      pb,&  ! 
                                      tk,&  ! temperature  
                                      theta,&! potential temperature 
                                      qv,&  ! water vapor 
                                      u,v   ! wind components 
! REAL,ALLOCATABLE,DIMENSION(:)    :: surf_values
  REAL    :: stand_lon,true_lat1,true_lat2,cen_lat,cen_lon,delta_x, &
             delta_y,delta_t
  REAL    :: min_lon,max_lon,range_lon,min_lat,max_lat,range_lat
  INTEGER :: map_projection

  INTEGER :: xlat_flg,xlon_flg,z_flg,p_flg,ph_flg,phb_flg,pp_flg, &
             pb_flg,tk_flg,theta_flg,qv_flg,parms_flg,u_flg,v_flg

  PRIVATE :: xlat,xlon,z,p,ph,phb,pp,pb,tk,theta,qv,u,v,          &
             xlat_flg,xlon_flg,z_flg,p_flg,ph_flg,phb_flg,pp_flg, &
             pb_flg,tk_flg,theta_flg,qv_flg,parms_flg,u_flg,v_flg
CONTAINS 
  SUBROUTINE init_module_wrf_getvar 
    xlat_flg = -1
    xlon_flg = -1 
    z_flg    = -1
    p_flg    = -1
    ph_flg   = -1 
    phb_flg  = -1
    pp_flg   = -1
    pb_flg   = -1 
    tk_flg   = -1
    theta_flg= -1
    qv_flg   = -1
    parms_flg= -1
    u_flg    = -1
    v_flg    = -1
  END SUBROUTINE init_module_wrf_getvar 

  SUBROUTINE getflg(var,flg)
    character(len=*), intent(in) :: var
    integer,      intent(out):: flg
    
    flg = LEN_TRIM(var)
    SELECT CASE (var(1:flg)) 
    CASE ('xlat_flg')
      flg = xlat_flg 
    CASE ('xlon_flg')
      flg = xlon_flg 
    CASE ('z_flg')
      flg = z_flg 
    CASE ('p_flg')
      flg = p_flg 
    CASE ('ph_flg')
      flg = ph_flg 
    CASE ('phb_flg')
      flg = phb_flg 
    CASE ('pp_flg')
      flg = pp_flg 
    CASE ('pb_flg')
      flg = pb_flg 
    CASE ('tk_flg')
      flg = tk_flg 
    CASE ('theta_flg')
      flg = theta_flg 
    CASE ('qv_flg')
      flg = qv_flg 
    CASE ('parms_flg')
      flg = parms_flg 
    CASE ('u_flg')
      flg = u_flg 
    CASE ('v_flg')
      flg = v_flg 
    CASE DEFAULT 
      write(6,'(2a)') ' The requested flag not exists -- ',var(1:flg)
      flg = -1
    END SELECT 
  END SUBROUTINE getflg 

  SUBROUTINE nullify_module_wrf_getvar
    if (allocated(r_buf)) deallocate(r_buf)
    if (allocated(xlon)) deallocate(xlon)
    if (allocated(xlat)) deallocate(xlat)
    if (allocated(p)) deallocate(p)
    if (allocated(pp)) deallocate(pp)
    if (allocated(pb)) deallocate(pb)
    if (allocated(z)) deallocate(z)
    if (allocated(ph)) deallocate(ph)
    if (allocated(phb)) deallocate(phb)
    if (allocated(qv)) deallocate(qv)
    if (allocated(theta)) deallocate(theta)
    if (allocated(tk)) deallocate(tk)
  END SUBROUTINE nullify_module_wrf_getvar

SUBROUTINE wrf_user_getvar(cdfid,varname,time,nx_out,ny_out,nz_out)
IMPLICIT NONE 

real               :: eps

INTEGER,          INTENT(in)   :: cdfid, time 
CHARACTER(len=*) ,INTENT(in)   :: varname
INTEGER,          INTENT(out)  :: nx_out, ny_out, nz_out
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: var,w
REAL, ALLOCATABLE, DIMENSION(:,:)   :: tmp1,t_surf,diff,alpha 
INTEGER, ALLOCATABLE, DIMENSION(:,:):: itmp1 

CHARACTER*12 :: var_to_get
REAL ::  cone,true_lat1,true_lat2,stand_lon,cen_lat,cen_lon
REAL ::  tmp_val, tmp_alpha
INTEGER :: i,j,k,map_projection,len_var,dims(4),ndims

  if( (varname .eq. "umet") .or. (varname .eq. "vmet") .or. & 
      (varname(1:4) .eq. "wind") ) then

    if (varname .eq. "wind10") then
      !the "V10" variable is staggered in north-south direction 
      var_to_get = "V10"
      len_var=len_trim(var_to_get)
      if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
      call get_wrf_v10(cdfid,time)
  
      !the "U10" variable is staggered in west-east direction 
      var_to_get = "U10"
      len_var=len_trim(var_to_get)
      if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
      call get_wrf_u10(cdfid,time)
  
      nx_out = nx
      ny_out = ny
      nz_out = 1
    else 
      !the "V" variable is staggered in north-south direction 
      var_to_get = "V"
      len_var=len_trim(var_to_get)
      if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
      call get_wrf_v(cdfid,time)
  
      !the "U" variable is staggered in west-east direction 
      var_to_get = "U"
      len_var=len_trim(var_to_get)
      if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
      call get_wrf_u(cdfid,time)
  
      nx_out = nx
      ny_out = ny
      nz_out = nz
    endif 
!!!
!!!  get projection's information and the value of the 'umet' or 'vmet'. 
!!!  
    call get_gl_att_int_cdf( cdfid, 'MAP_PROJ', map_projection, debug )

    if(map_projection .eq. 0) then  ! no projection
      if (allocated(r_buf)) deallocate(r_buf)
      allocate( r_buf(nx_out,ny_out,nz_out) )
      if(varname .eq. "vmet") then
        do i = 1,nx_out
          r_buf(i,:,:) = v(i,:,:)
        end do 
      else  ! must want umet
        do j = 1,ny_out
          r_buf(:,j,:) = u(:,j,:)
        enddo 
      end if

      return 
    end if

    call get_gl_att_real_cdf( cdfid, 'TRUELAT1', true_lat1, debug )
    call get_gl_att_real_cdf( cdfid, 'TRUELAT2', true_lat2, debug )

    cone = 1.
    if(map_projection .eq. 1) then    ! Lambert Conformal mapping
      if( (abs(true_lat1 - true_lat2) .gt. 0.1) .and. & 
        (abs(true_lat2 - 90. )      .gt. 0.1)       ) then
        cone = 10**(cos(true_lat1*radians_per_degree)) &
               -10**(cos(true_lat2*radians_per_degree))

!   ??? there may be errors in the function 'tan' of the 
!       following sentence 
!
        cone = cone/(10**(tan((45. -abs(true_lat1*0.5))*radians_per_degree)) &
               - 10**(tan((45. -abs(true_lat2*0.5))*radians_per_degree))   )
      else
        cone = sin(abs(true_lat1)*radians_per_degree)
      end if
    end if
    if(map_projection .eq. 2) then      ! polar steraographic
      cone = 1.
    end if
    if(map_projection .eq. 3) then      ! Mercator
      cone = 0.
    end if

!!!
!!! get latitude and longitude values of each grads. 
!!! 
  
    call get_wrf_xlat(cdfid,1) 
    call get_wrf_xlon(cdfid,1)

!   modified by Deng Liantang 2004.07.02
!   由于WRFOUT的输出结果中，嵌套网格的中心点位置没有给出，还是粗
!   网格的信息，所以在获取中心点位置和经纬都信息时，采用以下方法。
!   改正了用于计算离标准经度之间的经度偏差的标准经度值错误，正确的
!   值应该是stand_lon，而不是cen_lon. Deng Liantang 2006.07.10 
 
!   call get_gl_att_real_cdf( cdfid, 'CEN_LAT', cen_lat, debug )
!   call get_gl_att_real_cdf( cdfid, 'CEN_LON', cen_lon, debug )
!   cen_lat = (xlat(INT(nx/2.+0.5),INT(ny/2.+0.5)) &
!             +xlat(INT(nx/2.+1.0),INT(ny/2.+1.0)))/2.
!   cen_lon = (xlon(INT(nx/2.+0.5),INT(ny/2.+0.5)) &
!             +xlon(INT(nx/2.+1.0),INT(ny/2.+1.0)))/2.

      
    call get_gl_att_real_cdf( cdfid, 'STAND_LON', stand_lon, debug )
    allocate( diff(nx_out,ny_out) )
    diff = xlon - stand_lon

    do j = 1,ny_out 
      do i = 1,nx_out 
        if(diff(i,j) .gt. 180.) then
          diff(i,j) = diff(i,j) - 360.
        end if
        if(diff(i,j) .lt. -180.) then
          diff(i,j) = diff(i,j) + 360.
        end if
      end do
    end do

    allocate( alpha(nx_out,ny_out) )
!   alpha = diff * cone * radians_per_degree *sign(1.,latitude)
    do j = 1, ny_out
      do i = 1, nx_out
        if(xlat(i,j) .lt. 0.) then
          alpha(i,j) = - diff(i,j) * cone * radians_per_degree
        else
          alpha(i,j) = diff(i,j) * cone * radians_per_degree
        end if
      end do
    end do

    allocate( var(nx,ny,nz) )
    var=0.0
    if (allocated(r_buf)) deallocate(r_buf)
    allocate( r_buf(nx,ny,nz) )
    r_buf=0.0
      do k=1,nz_out
        do j=1,ny_out
          do i=1,nx_out
            r_buf(i,j,k) = v(i,j,k)*sin(alpha(i,j))+u(i,j,k)*cos(alpha(i,j))
            var(i,j,k) = v(i,j,k)*cos(alpha(i,j))-u(i,j,k)*sin(alpha(i,j))
          end do
        end do
      end do

    u=r_buf ; v=var 
    if(varname .eq. "vmet") then
      r_buf=v
    end if

    deallocate(var)
    deallocate(diff)
    deallocate(alpha)
    return 

  end if

  if( varname .eq. "ua" ) then

!   the "U" variable is staggered in west-east direction

    var_to_get = "U"
    len_var=len_trim(var_to_get)
    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
    call get_wrf_u(cdfid,time)
    
    nx_out = nx
    ny_out = ny
    nz_out = nz
    
!   copy the value of 3d variable u into output 3d variable r_buf .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = u
    return 

  end if

  if( varname .eq. "u" ) then
    var_to_get = "U"
    len_var=len_trim(var_to_get)
    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
    call get_wrf_u(cdfid,time)

    nx_out = nx
    ny_out = ny
    nz_out = nz

!   copy the value of 3d variable u into output 3d variable r_buf .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = u 
    return

  end if

  if( varname .eq. "va" ) then

!   the "V" variable is staggered in north-south direction

    var_to_get = "V"
    len_var=len_trim(var_to_get)
    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
    call get_wrf_v(cdfid,time)

!   copy the value of 3d variable v into output 3d variable r_buf .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx,ny,nz))
    r_buf = v 

    nx_out = nx
    ny_out = ny 
    nz_out = nz
    return

  end if

  if( varname .eq. "v" ) then
    var_to_get = "V"
    len_var=len_trim(var_to_get)
    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)
    call get_wrf_v(cdfid,time)

    nx_out = nx
    ny_out = ny
    nz_out = nz

!   copy the value of 3d variable v into output 3d variable r_buf .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = v

    return
  end if

  if( varname .eq. "wa" ) then

!   the "W" is staggered in vertical direction .

    var_to_get = "W"
    len_var=len_trim(var_to_get)

    dims(1) = nx; dims(2) = ny; dims(3) = nz_stg
    allocate( var(dims(1), dims(2), dims(3)) )

    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)

    call get_var_3d_real_cdf( cdfid, var_to_get(1:len_var), var, &
                              dims(1), dims(2), dims(3), time, debug)
    allocate( w(dims(1), dims(2), dims(3)-1) )
    do k=1,dims(3)-1
      w(:,:,k) = 0.5*(var(:,:,k) + var(:,:,k+1))
    enddo

    deallocate(var)

    nx_out = dims(1)
    ny_out = dims(2)
    nz_out = dims(3) - 1

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = w

    deallocate(w)
    return

  end if

  if( varname .eq. "w" ) then
    var_to_get = "W"
    len_var=len_trim(var_to_get)

    dims(1) = nx; dims(2) = ny; dims(3) = nz_stg
    allocate( var(dims(1), dims(2), dims(3)) )

    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)

    call get_var_3d_real_cdf( cdfid, var_to_get(1:len_var), var, &
                              dims(1), dims(2), dims(3), time, debug)

    nx_out = dims(1) 
    ny_out = dims(2)
    nz_out = dims(3)

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = var

    deallocate(var)
    return

  end if

  if( varname .eq. "om" ) then
! stagger in Z direction 
    var_to_get = "WW"
    len_var=len_trim(var_to_get)
    dims(1) = nx; dims(2) = ny; dims(3) = nz_stg
    allocate( var(dims(1), dims(2), dims(3)) )

    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)

    call get_var_3d_real_cdf( cdfid, var_to_get(1:len_var), var, &
                              dims(1), dims(2), dims(3), time, debug)

    do k = 1, dims(3) - 1
      var(:,:,k) = 0.5*(var(:,:,k)+var(:,:,k+1))
    enddo 
    nx_out = dims(1) 
    ny_out = dims(2)
    nz_out = dims(3)-1

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    do k = 1, dims(3) - 1
      r_buf(:,:,k) =var(:,:,k)
    enddo 

    deallocate(var)
    return

  end if

  if( varname .eq. "th" ) then
    call get_wrf_theta(cdfid,time)

    dims(1) = nx; dims(2) = ny; dims(3) = nz
    nx_out = dims(1) 
    ny_out = dims(2)
    nz_out = dims(3)

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = theta
     
    return 

  end if

  if( varname .eq. "p" ) then
    dims(1) = nx; dims(2) = ny; dims(3) = nz
    call get_wrf_p(cdfid,time)
    nx_out = dims(1) 
    ny_out = dims(2)
    nz_out = dims(3)

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = 0.01*p
 
    return 

  end if

  if( varname .eq. "Z" ) then
    dims(1) = nx; dims(2) = ny; dims(3) = nz
    call get_wrf_z(cdfid,time)

    nx_out = dims(1)
    ny_out = dims(2)
    nz_out = dims(3) 

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = z

    return

  end if

  if( varname .eq. "tc" ) then

    dims(1) = nx; dims(2) = ny; dims(3) = nz
    call get_wrf_tk(cdfid,time)

    nx_out = dims(1)
    ny_out = dims(2)
    nz_out = dims(3)

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = tk - t0_abs

    return

  end if

  if( varname .eq. "TSK" ) then

    var_to_get = "TSK"
    len_var=len_trim(var_to_get)

    dims(1) = nx; dims(2) = ny; dims(3) = 1
    if (allocated(r_buf)) deallocate(r_buf)
    allocate( r_buf(dims(1), dims(2), dims(3)) )

    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)

    call get_var_2d_real_cdf( cdfid, var_to_get(1:len_var), r_buf, &
                              dims(1), dims(2), time, debug)
    r_buf = r_buf - t0_abs

    nx_out = dims(1)
    ny_out = dims(2)
    nz_out = dims(3)
    return

  end if 

  if( varname .eq. "T2" ) then

    var_to_get = "T2"
    len_var=len_trim(var_to_get)

    dims(1) = nx; dims(2) = ny; dims(3) = 1
    if (allocated(r_buf)) deallocate(r_buf)
    allocate( r_buf(dims(1), dims(2), dims(3)) )

    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)

    call get_var_2d_real_cdf( cdfid, var_to_get(1:len_var), r_buf, &
                              dims(1), dims(2), time, debug)
    r_buf = r_buf - t0_abs

    nx_out = dims(1)
    ny_out = dims(2)
    nz_out = dims(3)
    return

  end if 

  if( varname .eq. "rain" ) then

    var_to_get = "RAINC"
    len_var=len_trim(var_to_get)

    dims(1) = nx; dims(2) = ny; dims(3) = 1

    if (allocated(r_buf)) deallocate(r_buf)
    allocate( r_buf(dims(1), dims(2), dims(3)) )

    call get_var_2d_real_cdf( cdfid, var_to_get(1:len_var), r_buf, &
                              dims(1), dims(2), time, debug)

    var_to_get = "RAINNC"
    len_var=len_trim(var_to_get)

    if (allocated(tmp1)) deallocate(tmp1)
    allocate( tmp1(dims(1), dims(2) ) )

    if(debug .GT. 0) write(6,*)' getting data for ',var_to_get(1:len_var)

    call get_var_2d_real_cdf( cdfid, var_to_get(1:len_var), tmp1, &
                              dims(1), dims(2), time, debug)

!   r_buf(:,:,1) = r_buf(:,:,1) + tmp1
    tmp1 = r_buf(:,:,1) + tmp1 

    nx_out = dims(1)
    ny_out = dims(2)
    nz_out = dims(3)

    ! 为了处理降水量而添加了一个r_save的量，保存前一时次的降水量 2005.06.21
    ! 有逻辑错误，当读完其他变量后再读rain时，保存的前一时次量就不是降水，
    ! 量了,而是其他量。
    ! 20060530 已经修正了该逻辑错误。  
    ! r_save 是前一时次累积总量，r_buf是逐时间帧降水量(如6小时输出一次，则是6小时累积量。)

    if (.NOT. allocated(r_save)) then 
      r_buf(:,:,1) = tmp1 
      allocate(r_save(dims(1), dims(2)))
      r_save = 0.0 
    else 
      r_buf(:,:,1) = tmp1 - r_save
      r_save = tmp1 
    endif 

    deallocate(tmp1)

    return

  end if 

  if( varname .eq. "td" ) then

    ! compute qv, p, then td
    dims(1) = nx; dims(2) = ny; dims(3) = nz
    call get_wrf_qv(cdfid,time)
    where (qv < 0.0) qv = 0.0 
    call get_wrf_p(cdfid,time)

    allocate(var(dims(1),dims(2),dims(3)) ) 

    var = qv * 0.01*p / ( 0.622+qv) ! vapor pressure 

    where (var < 0.01) var = 0.01       ! avoid problems near zero 

!   var = (243.5/( (17.67/log(var/6.112)) - 1.0)) ! Bolton's 
                                                  ! approximation
    var = (243.5*log(var)-440.8)/(19.48-log(var))

    nx_out = dims(1) 
    ny_out = dims(2)
    nz_out = dims(3)

    ! copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = var
    deallocate(var)
    return

  end if

  if( varname .eq. "td2" ) then

    nx_out = nx
    ny_out = ny
    nz_out = 1

    if (.NOT. allocated(r_buf)) allocate(r_buf(nx_out,ny_out,nz_out))
    call compute_td2(cdfid,time,r_buf)

    return

  end if

  if( varname .eq. "slvl" ) then

! compute tk, p, qv, z, then slvl  
    dims(1) = nx; dims(2) = ny; dims(3) = nz
    call get_wrf_tk(cdfid,time)
    call get_wrf_qv(cdfid,time)
    where (qv < 0.0) qv = 0.0
    call get_wrf_z(cdfid,time)
    call get_wrf_p(cdfid,time)

    allocate( var(dims(1),dims(2),1) )
    allocate( t_surf(dims(1),dims(2)) )
    allocate( tmp1(dims(1),dims(2)) )
    allocate( itmp1(dims(1),dims(2)) )

    call compute_seaprs( cdfid,time, var, t_surf,tmp1, itmp1 )

    var = 0.01*var

    nx_out = dims(1) 
    ny_out = dims(2)
    nz_out = 1

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = var 
    
    deallocate(var)
    deallocate(t_surf)
    deallocate(tmp1) 
    deallocate(itmp1)

    return

  end if

  if( varname .eq. "rh" ) then

! compute qv, p, tk, then rh

    call get_wrf_tk(cdfid,time)
    call get_wrf_qv(cdfid,time)
    where (qv < 0.0) qv = 0.0
    call get_wrf_p(cdfid,time)

    allocate( var(nx,ny,nz) )
    call compute_rh(cdfid,time,var)

    nx_out = nx
    ny_out = ny
    nz_out = nz

!   copy the value of 3d variable "var" into output 3d variable "data" .
    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = var

    deallocate(var) 
    return

  end if

  if( varname .eq. "rh2" ) then

    nx_out = nx
    ny_out = ny
    nz_out = 1 

    if (.NOT. allocated(r_buf)) allocate(r_buf(nx_out,ny_out,nz_out))
    call compute_rh2(cdfid,time,r_buf)

    return

  end if

!  end of diagnostic varname list 
  
!  We want other fields already in this file. Check variable 
!  dimensionality and pull proper time out of file.

    var_to_get = varname
    len_var=len_trim(var_to_get)

    if(debug .GT. 0) write(6,*)' getting dims for ',var_to_get(1:len_var)

    call get_dims_cdf( cdfid, var_to_get(1:len_var), &
                       dims, ndims, debug) 

    SELECT CASE (ndims ) 
    CASE (4)
      allocate( var(dims(1), dims(2), dims(3)) )

      if(debug .GT. 0) write(6,*)' 4 dims: getting data for ',var_to_get(1:len_var)

      call get_var_3d_real_cdf( cdfid, var_to_get(1:len_var), var, &
                                dims(1), dims(2), dims(3), time, debug)

      nx_out = dims(1)
      ny_out = dims(2)
      nz_out = dims(3)

!   copy the value of 3d variable "var" into output 3d variable "data" .
      if (allocated(r_buf)) deallocate(r_buf)
      allocate(r_buf(nx_out,ny_out,nz_out))
      r_buf = var

      deallocate(var)

    CASE (3) 

      allocate( var(dims(1), dims(2), 1) )

      if(debug .GT. 0) write(6,*)' 3 dims:getting data for ',var_to_get(1:len_var)

      call get_var_2d_real_cdf( cdfid, var_to_get(1:len_var), var, &
                                dims(1), dims(2), time, debug)

      nx_out = dims(1)
      ny_out = dims(2)
      nz_out = 1

!   copy the value of 3d variable "var" into output 3d variable "data" .
      if (allocated(r_buf)) deallocate(r_buf)
      allocate(r_buf(nx_out,ny_out,nz_out))
      r_buf = var
 
      deallocate(var)
    CASE (2) 

      allocate( var(dims(1), 1, 1) )

      if(debug .GT. 0) write(6,*)' 2 dims:getting data for ',var_to_get(1:len_var)

      call get_var_1d_int_cdf( cdfid, var_to_get(1:len_var), var, &
                                dims(1), time, debug)

      nx_out = dims(1)
      ny_out = 1
      nz_out = 1

!   copy the value of 3d variable "var" into output 3d variable "data" .
      if (allocated(r_buf)) deallocate(r_buf)
      allocate(r_buf(nx_out,ny_out,nz_out))
      r_buf = var
 
      deallocate(var)

    CASE (1) 

    ! now, we don't deal with the 1D variables
      nx_out = 1; ny_out = 1; nz_out = 1 
      if (allocated(r_buf)) deallocate(r_buf)
      allocate(r_buf(nx_out,ny_out,nz_out))
      r_buf = miss_val

    CASE DEFAULT 
      write(6,*) " The dimensions of variable ",varname(1:len_var), &
                 "(",ndims,") is wrong!" 
      stop " Wrong dimensions" 
    END SELECT 
    
    return 

END SUBROUTINE wrf_user_getvar

SUBROUTINE get_wrf_model_parms(cdfid) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid

  if (parms_flg .EQ. 1) return
  if(debug .GE. 100) write(6,*)' get_wrf_model_parms: getting model global parameters'

! 加密输出时，有时只有区域的常规大小信息，没有跳格点的大小信息，故采用
! 常规区域大小来推算跳格点的区域大小信息。

  !call get_gl_dims_len_cdf(cdfid,"west_east_stag",nx_stg,debug)
  !call get_gl_dims_len_cdf(cdfid,"south_north_stag",ny_stg,debug)
  !call get_gl_dims_len_cdf(cdfid,"bottom_top_stag",nz_stg,debug)

  call get_gl_dims_len_cdf(cdfid,"west_east",nx,debug)
  call get_gl_dims_len_cdf(cdfid,"south_north",ny,debug)
  call get_gl_dims_len_cdf(cdfid,"bottom_top",nz,debug)
  nx_stg = nx + 1; ny_stg = ny + 1; nz_stg = nz + 1   

  call get_gl_dims_len_cdf(cdfid,"soil_layers_stag",soil_stg,debug)
  call get_gl_att_int_cdf( cdfid, 'MAP_PROJ', map_projection, debug )
  call get_gl_att_real_cdf( cdfid, 'TRUELAT1', true_lat1, debug )
  call get_gl_att_real_cdf( cdfid, 'TRUELAT2', true_lat2, debug )
  call get_gl_att_real_cdf( cdfid, 'STAND_LON',stand_lon, debug )

  call get_gl_att_real_cdf( cdfid, 'CEN_LAT', cen_lat, debug )
  call get_gl_att_real_cdf( cdfid, 'CEN_LON', cen_lon, debug )

  call get_gl_att_real_cdf( cdfid, 'DX', delta_x, debug )
  call get_gl_att_real_cdf( cdfid, 'DY', delta_y, debug )
  call get_gl_att_real_cdf( cdfid, 'DT', delta_t, debug )

  !; stand_lon = cen_lon   ! 在嵌套网格中，stand_lon 可能与cen_lon不完全一样 
    
! modified by Deng Liantang 2004.07.02
! 由于WRFOUT的输出结果中，嵌套网格的中心点位置没有给出，还是粗
! 网格的信息，所以在获取中心点位置和经纬都信息时，采用以下方法。
!   call get_wrf_xlat(cdfid,1)
!   call get_wrf_xlon(cdfid,1)
!   cen_lat = (xlat(INT(nx/2.+0.5),INT(ny/2.+0.5)) &
!             +xlat(INT(nx/2.+1.0),INT(ny/2.+1.0)))/2.
!   cen_lon = (xlon(INT(nx/2.+0.5),INT(ny/2.+0.5)) &
!             +xlon(INT(nx/2.+1.0),INT(ny/2.+1.0)))/2 

  parms_flg = 1

END SUBROUTINE get_wrf_model_parms

SUBROUTINE get_wrf_pp(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if (pp_flg .EQ. time) return
  if (.NOT. allocated(pp)) allocate(pp(nx,ny,nz))

  if(debug .GE. 100) write(6,*)' get_wrf_pp: getting data for P'
  call get_var_3d_real_cdf( cdfid, "P" , pp, nx, ny, nz, time, debug)
  
  pp_flg = time

END SUBROUTINE get_wrf_pp

SUBROUTINE get_wrf_pb(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if (pb_flg .EQ. time) return
  if (.NOT. allocated(pb)) allocate(pb(nx,ny,nz))

  if(debug .GE. 100) write(6,*)' get_wrf_pb: getting data for PB'
  call get_var_3d_real_cdf( cdfid, "PB" , pb, nx, ny, nz, time, debug)
  
  pb_flg = time

END SUBROUTINE get_wrf_pb 

SUBROUTINE get_wrf_ph(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 
  
  if (ph_flg .EQ. time) return
  if (.NOT. allocated(ph)) allocate(ph(nx,ny,nz_stg))

  if(debug .GE. 100) write(6,*)' get_wrf_ph: getting data for PH'
  call get_var_3d_real_cdf( cdfid, "PH" , ph, nx, ny, nz_stg, time, debug)
  
  ph_flg = time

END SUBROUTINE get_wrf_ph 

SUBROUTINE get_wrf_phb(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if (phb_flg .EQ. time) return
  if (.NOT. allocated(phb)) allocate(phb(nx,ny,nz_stg))

  if(debug .GE. 100) write(6,*)' get_wrf_phb: getting data for PHB'
  call get_var_3d_real_cdf( cdfid, "PHB" , phb, nx, ny, nz_stg, time, debug)
  
  phb_flg = time

END SUBROUTINE get_wrf_phb 

SUBROUTINE get_wrf_qv(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if (qv_flg .EQ. time) return
  if (.NOT. allocated(qv)) allocate(qv(nx,ny,nz))

  if(debug .GE. 100) write(6,*)' get_wrf_qv: getting data for QVAPOR' 
  call get_var_3d_real_cdf( cdfid, "QVAPOR" , qv, nx, ny, nz, time, debug)
  
  qv_flg = time

END SUBROUTINE get_wrf_qv 

SUBROUTINE get_wrf_theta(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if (theta_flg .EQ. time) return
  if (.NOT. allocated(theta)) allocate(theta(nx,ny,nz))

  if(debug .GE. 100) write(6,*)' get_wrf_theta: getting data for T ' 
  call get_var_3d_real_cdf( cdfid, "T" , theta, nx, ny, nz, time, debug)
  theta = theta + 300.
  
  theta_flg = time

END SUBROUTINE get_wrf_theta 

SUBROUTINE get_wrf_xlon(cdfid, time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if (xlon_flg .EQ. time) return
  if (.NOT. allocated(xlon)) allocate(xlon(nx,ny))

  if(debug .GE. 100) write(6,*)' get_wrf_xlon: getting data for XLONG ' 
  call get_var_2d_real_cdf( cdfid, "XLONG", xlon, nx, ny, time, debug)
  min_lon=MINVAL(xlon); max_lon=MAXVAL(xlon); range_lon=max_lon-min_lon
  if(debug .GE. 100) write(6,*)' min_lon/max_lon/rangle_lon: ',min_lon,max_lon,range_lon  

  xlon_flg = time

END SUBROUTINE get_wrf_xlon 

SUBROUTINE copy_wrfvar(varname,nx_out,ny_out,nz_out) 
  IMPLICIT NONE
  CHARACTER(len=*), intent(in)  :: varname
  integer,          intent(out) :: nx_out,ny_out,nz_out

  INTEGER  :: len_var,dims(3)
 
  if (allocated(r_buf)) deallocate(r_buf)
  len_var = len_trim(varname)

  SELECT CASE (varname(1:len_var)) 
  CASE ("Z","z")
    dims=shape(z)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = z
  CASE ("P","p")
    dims=shape(p)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = p
  CASE ("tk","TK")
    dims=shape(tk)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = tk
  CASE ("qv","QV")
    dims=shape(qv)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = qv
  CASE ("theta","THETA")
    dims=shape(theta)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = theta
  CASE ("pp","PP")
    dims=shape(pp)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = pp
  CASE ("pb","PB")
    dims=shape(pb)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = pb
  CASE ("ph","PH")
    dims=shape(ph)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = ph
  CASE ("phb","PHB")
    dims=shape(phb)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = phb
  CASE ("u","U")
    dims=shape(u)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = u
  CASE ("v","V")
    dims=shape(v)
    nx_out = dims(1); ny_out = dims(2); nz_out = dims(3)
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf = v
  CASE ("XLAT","xlat")
    dims(1:2)=shape(xlat)
    nx_out = dims(1); ny_out = dims(2); nz_out = 1
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf(:,:,1) = xlat
  CASE ("XLON","xlon")
    dims(1:2)=shape(xlon)
    nx_out = dims(1); ny_out = dims(2); nz_out = 1
    allocate(r_buf(nx_out,ny_out,nz_out))
    r_buf(:,:,1) = xlon
  CASE DEFAULT 
    write(6,'(25("*"),"ERROR",30("*"))')
    write(6,*) " copy_wrfvar: can not copy the r_bufiable: ",varname(1:len_var)
    write(6,'(60("*"))')
  END SELECT 

  return 
END SUBROUTINE copy_wrfvar

SUBROUTINE get_wrf_xlat(cdfid,time) 
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

  if ( xlat_flg .EQ. time) return
  if (.NOT. allocated(xlat)) allocate(xlat(nx,ny))

  if(debug .GE. 100) write(6,*)' get_wrf_xlat: getting data for XLAT ' 
  call get_var_2d_real_cdf( cdfid, "XLAT", xlat, nx, ny, time, debug)
  min_lat=MINVAL(xlat); max_lat=MAXVAL(xlat); range_lat=max_lat-min_lat
  if(debug .GE. 100) write(6,*)' min_lat/max_lat/rangle_lat: ',min_lat,max_lat,range_lat  

  xlat_flg = time

END SUBROUTINE get_wrf_xlat 

SUBROUTINE get_wrf_p(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 
  INTEGER                :: i,j,k

    if (p_flg .EQ. time) return
    if (.NOT. allocated(p)) allocate(p(nx,ny,nz))
    
    if(debug .GE. 100) write(6,*)' get_wrf_p: getting data for P and PB ' 
    call get_wrf_pp(cdfid,time)
    call get_wrf_pb(cdfid,time)

    do k=1,nz; do j=1,ny; do i=1,nx 
      p(i,j,k) = pp(i,j,k) + pb(i,j,k) 
    enddo ; enddo; enddo 

    p_flg = time

END SUBROUTINE get_wrf_p

SUBROUTINE get_wrf_z(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 
  INTEGER                :: i,j,k
  REAL                   :: g1

    if (z_flg .EQ. time) return
    if (.NOT. allocated(z)) allocate(z(nx,ny,nz))
    
    if(debug .GE. 100) write(6,*)' get_wrf_z: getting data for PH and PHB ' 
!   the "PH" variable is staggered in vertical direction 
    call get_wrf_ph(cdfid,time)

!   the "PHB" variable is staggered in vertical direction. 
    call get_wrf_phb(cdfid,time)

    g1=0.5/g
    do k=1,nz; do j=1,ny; do i=1,nx 
      z(i,j,k) = g1*(ph(i,j,k)+ph(i,j,k+1) + phb(i,j,k)+phb(i,j,k+1))
    enddo ; enddo; enddo 

    z_flg = time

END SUBROUTINE get_wrf_z

SUBROUTINE get_wrf_tk(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

    if (tk_flg .EQ. time) return
    if(debug .GE. 100) write(6,*)' get_wrf_tk: getting data for T, P, PB' 

! compute theta, p, and then tc

    call get_wrf_theta(cdfid,time)
    call get_wrf_p(cdfid,time)
    call compute_tk(cdfid,time)
    
    tk_flg = time

END SUBROUTINE get_wrf_tk

SUBROUTINE get_wrf_v(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 
  INTEGER :: i,j,k 

    if (v_flg .EQ. time) return

    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx,ny_stg,nz))
    if (.NOT. allocated(v)) allocate(v(nx,ny,nz))

! "V" is staggered in south-north direction 

    if(debug .GE. 100) write(6,*)' get_wrf_v: getting data for V' 
    call get_var_3d_real_cdf(cdfid,"V",r_buf,nx,ny_stg,nz,time,debug)

    do k=1,nz; do j=1,ny; do i=1,nx 
      v(i,j,k) = 0.5*(r_buf(i,j,k)+r_buf(i,j+1,k))
    enddo ; enddo; enddo 

    v_flg = time

END SUBROUTINE get_wrf_v

SUBROUTINE get_wrf_v10(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

    ! "V10" is not staggered in south-north direction 
    if (.NOT. allocated(v)) allocate(v(nx,ny,nz))
    v=0.0

    if(debug .GE. 100) write(6,*)' get_wrf_v10: getting data for V10' 
    call get_var_2d_real_cdf(cdfid,"V10",v(:,:,1),nx,ny,time,debug)

    v_flg = -1

END SUBROUTINE get_wrf_v10

SUBROUTINE get_wrf_u(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 
  INTEGER :: i,j,k 

    if (u_flg .EQ. time) return

    if (allocated(r_buf)) deallocate(r_buf)
    allocate(r_buf(nx_stg,ny,nz))
    if (.NOT. allocated(u)) allocate(u(nx,ny,nz))

! "U" is staggered in west-east direction 

    if(debug .GE. 100) write(6,*)' get_wrf_u: getting data for U' 
    call get_var_3d_real_cdf(cdfid,"U",r_buf,nx_stg,ny,nz,time,debug)

    do k=1,nz; do j=1,ny; do i=1,nx 
      u(i,j,k) = 0.5*(r_buf(i,j,k)+r_buf(i+1,j,k))
    enddo ; enddo; enddo 

    u_flg = time

END SUBROUTINE get_wrf_u

SUBROUTINE get_wrf_u10(cdfid,time)
  IMPLICIT NONE
  INTEGER,  INTENT(in)   :: cdfid, time 

    ! "U10" is not staggered in west-east direction 
    if (.NOT. allocated(u)) allocate(u(nx,ny,nz))
    u=0.0

    if(debug .GE. 100) write(6,*)' get_wrf_u10: getting data for U10' 
    call get_var_2d_real_cdf(cdfid,"U10",u(:,:,1),nx,ny,time,debug)

    u_flg = -1

END SUBROUTINE get_wrf_u10

SUBROUTINE get_u(u_out)
  IMPLICIT NONE
  INTEGER  :: i,j,k
  real,     INTENT(out)   :: u_out(nx,ny,nz)

  do k=1,nz; do j=1,ny; do i=1,nx 
    u_out(i,j,k) = u(i,j,k) 
  enddo ; enddo; enddo 
END SUBROUTINE get_u

SUBROUTINE get_v(v_out)
  IMPLICIT NONE
  real,     INTENT(out)   :: v_out(nx,ny,nz)
  INTEGER  :: i,j,k

  do k=1,nz; do j=1,ny; do i=1,nx 
    v_out(i,j,k) = v(i,j,k) 
  enddo ; enddo; enddo 
END SUBROUTINE get_v

       SUBROUTINE wrf_get_interp_3dvar(cdfid,varname,time, &
                                       interp,nlevel,layers,&
                                       nx_out,ny_out,nz_out)
       implicit none 

       integer,         intent(in) :: cdfid 
       integer,         intent(in) :: time  
       character(len=*),intent(in) :: varname
       integer,         intent(out):: nx_out, ny_out, nz_out
       integer,         intent(in) :: interp
       integer,         intent(in) :: nlevel
       real,            intent(in) :: layers(nlevel)

       real,allocatable :: var(:,:,:),var3d_tmp(:,:,:)
       real,allocatable :: var2d_out(:,:)
       real             :: levels(nlevel)
!       integer          :: ndims,dims(4)
       integer          :: k
       character*80     :: var_to_get
       character*1      :: p_str
       integer          :: len_var  
       integer          :: nlevel_out  

       levels = layers
       var_to_get = varname
       len_var=len_trim(var_to_get)

       ! get the value of 3D variable 
       call wrf_user_getvar(cdfid,var_to_get(1:len_var),time,  & 
                            nx_out,ny_out,nz_out)

       SELECT CASE (interp)
       CASE (0)
         p_str = ""
         nlevel_out = nz_out
       CASE (1)
         p_str = "p"
         nlevel_out = nlevel
         if (MAXVAL(levels) .LT. 3000.) then !maybe have the surface level 2001
           levels =  levels * 100.
         end if 
       CASE (2) 
         p_str = "Z"
         nlevel_out = nlevel
         if (MAXVAL(levels) .LT. 100.) then 
           levels =  levels * 1000.
         end if 
       CASE DEFAULT
         p_str = ""
         nlevel_out = nz_out
       END SELECT
       
       SELECT CASE (interp) 
       CASE (1,2)      ! need to interpolate variable 
         ! keep the obtained variable 
         allocate(var(nx_out,ny_out,nz_out)) 
         var = r_buf

         ! get the 3D variable p or z value used to interpolate
         if (p_str .eq. "p".or.p_str .eq. "P") then  
          call get_wrf_p(cdfid,time)
         elseif (p_str .eq. "z".or.p_str .eq. "Z") then
          call get_wrf_z(cdfid,time)
         endif 

         if (debug .GE. 100) write(6,*)" wrf_get_interp_3dvar: Interpolate variable ", &
                               var_to_get(1:len_var), " ..."

         ! un-staggered the dimensions size of var. 
         deallocate(r_buf)
         allocate(r_buf(nx_out,ny_out,nz_out))
         r_buf = var
         call xy_unstagger(nx_out,ny_out,nz_out)
         nx_out = nx; ny_out = ny
         call z_unstagger(nx_out,ny_out,nz_out)
         deallocate(var)

         ! interpolate the 3D variable in z direction 
         call interp_3dvar(p_str,levels,nlevel)
         nz_out = nlevel

       CASE (0)        ! not need to interpolate, just skip

       CASE DEFAULT ! CASE (interp)
         ierr = 1 
         write(6,*) ' Now, couldnot deal with this case: interp = ',interp 
         return 
       END SELECT 
       return

       END SUBROUTINE wrf_get_interp_3dvar

!  un-stagger the variable (store in r_buf) in model size (nx,ny,nz)
      SUBROUTINE z_unstagger( nx_in,ny_in,nz_in )
      implicit none

      integer,  intent(in) :: nx_in,ny_in,nz_in

      real                 :: newvar(nx_in,ny_in,nz)
      integer              :: i,j,k
      integer              :: n_count
      logical              :: staggered 

      ierr = 0 
      eps = ABS(epsilon * miss_val)
      staggered = .false.

      if (nz_in.LT.nz) then 
         if (debug .GT. 0) write(6,*) " ERROR z_unstagger: dimension mismatch" 
         if (debug .GT. 0) write(6,*) "       nz: ",nz,"nz_in: ",nz_in 
         ierr = 1
      endif 

      newvar = 0
      if(nz_in .gt. nz) then  ! w (z) staggering
        do k=1,nz
          do j=1,ny_in
            do i=1,nx_in
              n_count = 0
              if (ABS(r_buf(i,j,k)-miss_val) .GT. eps) then 
                newvar(i,j,k) = r_buf(i,j,k)
                n_count = n_count + 1
              end if 
              if (ABS(r_buf(i,j,k+1)-miss_val) .GT. eps) then 
                newvar(i,j,k) = newvar(i,j,k) + r_buf(i,j,k+1)
                n_count = n_count + 1
              end if 
              if (n_count .EQ. 2) then 
                newvar(i,j,k) = newvar(i,j,k) * 0.5 
              elseif (n_count .EQ. 0) then  
                newvar(i,j,k) = miss_val
              end if 
            enddo
          enddo
        enddo
        staggered = .true. 
      end if 
      if (staggered) then      ! un-staggered 
        deallocate(r_buf) 
        allocate(r_buf(nx_in,ny_in,nz))
        r_buf = newvar
      endif

      return
      END SUBROUTINE z_unstagger 

      SUBROUTINE  xy_unstagger(nx_in, ny_in, nz_in)
      implicit none 

      integer,      intent(in) :: nx_in, ny_in, nz_in

      real             :: newvar(nx,ny,nz_in) 
      integer          :: n_count
      integer          :: i,j,k
      logical          :: staggered 
      
      ierr = 0
      eps = ABS(epsilon * miss_val)
      staggered = .false.

      if (nx_in.LT.nx .or. ny_in.LT.ny) then 
         if(debug .GT. 0) write(6,*)" ERROR xy_unstagger: dimension mismatch" 
         if(debug .GT. 0) write(6,*)"  nx/ny: ",nx,ny,"nx_in/ny_in: ",nx_in,ny_in 
         ierr = 1
      endif 

      newvar = 0
      if (nx_in .GT. nx) then 
        do k = 1,nz_in 
          do j=1,ny
            do i=1,nx
              n_count = 0
              if (ABS(r_buf(i,j,k)-miss_val) .GT. eps) then
                newvar(i,j,k) = r_buf(i,j,k)
                n_count = n_count + 1
              end if
              if (ABS(r_buf(i+1,j,k)-miss_val) .GT. eps) then
                newvar(i,j,k) = newvar(i,j,k)+r_buf(i+1,j,k)
                n_count = n_count + 1
              end if
              if (n_count .EQ. 2) then
                newvar(i,j,k) = newvar(i,j,k) * 0.5
              elseif (n_count .EQ. 0) then
                newvar(i,j,k) = miss_val
              end if
            enddo
          enddo
        enddo
        if(debug .GE. 100) write(6,*) " unstaggered in x direction: ",nx_in,"->",nx
        staggered = .true. 
      end if 
      if (ny_in .GT. ny) then 
        do k=1,nz_in
          do j=1,ny
            do i=1,nx
              n_count = 0
              if (ABS(r_buf(i,j,k)-miss_val) .GT. eps) then
                newvar(i,j,k) = r_buf(i,j,k)
                n_count = n_count + 1
              end if
              if (ABS(r_buf(i,j+1,k)-miss_val) .GT. eps) then
                newvar(i,j,k) = newvar(i,j,k)+r_buf(i,j+1,k)
                n_count = n_count + 1
              end if
              if (n_count .EQ. 2) then
                newvar(i,j,k) = newvar(i,j,k) * 0.5
              elseif (n_count .EQ. 0) then
                newvar(i,j,k) = miss_val
              end if
            enddo
          enddo
        end do
        if (debug .GE. 100) write(6,*) " unstaggere in y direction: ", ny_in, "->", ny 
        staggered = .true.
      end if 
      if (staggered) then 
        deallocate(r_buf) 
        allocate(r_buf(nx,ny,nz_in))
        r_buf = newvar
      end if 
      END SUBROUTINE xy_unstagger 

      SUBROUTINE interp_3dvar(p_str, levels, nlevel)

      implicit none
      
      character(len=*),intent(in) :: p_str
      integer,         intent(in) :: nlevel
      real,            intent(in) :: levels(nlevel)

      real                :: var(nx,ny,nlevel),surf(nx,ny,nz)
      integer             :: i,j,k, kp, ip, im
      logical             :: interp
      real                :: height, w1, w2


! does vertical coordinate increase of decrease with increasing k?
! set offset appropriately

      ip = 0
      im = 1
      if (p_str .EQ. "p" .or. p_str .EQ. "P") then 
        surf = p 
      elseif (p_str .EQ. "z" .or. p_str .EQ. "Z") then  
        surf = z 
      else 
        ierr = 1
        if (debug .GE. 100) write(6,*) " ERROR interp_3dvar: Can not ", &
                              "deal with this surfaces: ",p_str
      endif 

      if (surf(1,1,1) .gt. surf(1,1,nz)) then
        ip = 1
        im = 0
      endif

      eps = ABS(miss_val*epsilon)
      DO k = 1, nlevel 
        height = levels(k) 
        DO j=1,ny
          DO i=1,nx
            interp = .false.
            kp = nz
            DO WHILE ( (.not. interp) .and. (kp .ge. 2) )

              IF(   ((surf(i,j,kp-im) .le. height) .and.        &
                    (surf(i,j,kp-ip) .gt. height))      )   THEN
               if ( (ABS(r_buf(i,j,kp-im)-miss_val).GT.eps).and. &
                         (ABS(r_buf(i,j,kp-ip)-miss_val).GT.eps) ) then 
                w2 = (height-surf(i,j,kp-im))/(surf(i,j,kp-ip)-surf(i,j,kp-im))
                w1 = 1.-w2
                var(i,j,k) = w1*r_buf(i,j,kp-im) + w2*r_buf(i,j,kp-ip)
                interp = .true.
               elseif (ABS(r_buf(i,j,kp-im)-miss_val).GT.eps) then 
                var(i,j,k) = r_buf(i,j,kp-im) 
                interp = .true.
               elseif (ABS(r_buf(i,j,kp-ip)-miss_val).GT.eps) then 
                var(i,j,k) = r_buf(i,j,kp-ip) 
                interp = .true.
               end if 
              END IF
              kp = kp-1

            ENDDO
               
            ! 赋成缺测值 
            if (.not. interp) then 
!yry              IF( (surf(i,j,1)*im + surf(i,j,nz)*ip) .gt. height ) THEN   
                ! 赋成模式最低层 
!yry                var(i,j,k) = r_buf(i,j,1)*im + r_buf(i,j,nz)*ip
!yry              ELSEIF( (surf(i,j,nz)*im + surf(i,j,1)*ip) .lt. height ) THEN   
                ! 赋成模式最高层 
!yry                var(i,j,k) = r_buf(i,j,nz)*im + r_buf(i,j,1)*ip
!yry              ELSE
                var(i,j,k) = miss_val
!yry              ENDIF 
            endif !if (.not. interp) 
              
          ENDDO 
        ENDDO
      ENDDO

      deallocate(r_buf)
      allocate(r_buf(nx,ny,nlevel)) 
      
      r_buf = var

      RETURN
      END SUBROUTINE interp_3dvar

SUBROUTINE compute_theta(density,base_theta)
implicit none

real, intent(in) ::    base_theta(nx,ny,nz)
real, intent(in) ::    density(nx,ny,nz)

integer i,j,k

! input theta is perturbation moist theta.

do k=1,nz
  do j=1,ny
    do i=1,nx
      theta(i,j,k) = (theta(i,j,k)+base_theta(i,j,k))/  &
                    (density(i,j,k)*(1.+1.61*qv(i,j,k)))
    enddo
  enddo
enddo

END SUBROUTINE compute_theta

      SUBROUTINE compute_tk(cdfid,time)
      implicit none
      integer , intent(in) :: cdfid, time

      real    pi
      integer i,j,k
      real p1000mb, r_d, cp
      parameter ( p1000mb = 100000., r_d = 287., cp = 7.*r_d/2. )

! get p and theta 
      call get_wrf_p(cdfid,time)
      call get_wrf_theta(cdfid,time)

      if (.NOT. allocated(tk)) allocate(tk(nx,ny,nz))

      do k=1,nz
        do j=1,ny
          do i=1,nx
            pi=(p(i,j,k)/p1000mb)**(r_d/cp)
            tk(i,j,k) = pi*theta(i,j,k)
          enddo
        enddo
      enddo

      return
      END SUBROUTINE compute_tk

!---------------------------------------------

! Bill,
! This routine assumes
!    index order is (i,j,k)
!    wrf staggering
!    units: pressure (Pa), temperature(K), height (m), mixing ratio (kg kg{-1})
!    availability of 3d p, t, and qv; 2d terrain; 1d half-level zeta string
!    output units of SLP are Pa, but you should divide that by 100 for the
!    weather weenies.
!    virtual effects are included
!
! I have done no testing on this routine.  Please feel free to give me a call
! on Monday  208 362-9747
! Dave

      SUBROUTINE compute_seaprs(cdfid,time,psml,tsml,t_surf,level)
      IMPLICIT NONE
!     Estimate sea level pressure.
      INTEGER,  intent(in)  :: cdfid, time

!     The output is the 2d sea level pressure.
      REAL,     intent(out) :: psml(nx,ny)
      INTEGER , intent(out) :: level(nx,ny)
      REAL,     intent(out) :: t_surf(nx,ny) , tsml(nx,ny)

!     Some required physical constants:

      REAL R, GAMMA
      PARAMETER (R=287.04, GAMMA=-0.0065)

!     Specific constants for assumptions made in this routine:

      REAL    TC, PCONST
      PARAMETER (TC=273.16+17.5, PCONST = 10000)
      LOGICAL ridiculous_mm5_test
      PARAMETER (ridiculous_mm5_test = .TRUE.)
!      PARAMETER (ridiculous_mm5_test = .false.)

!     Local variables:

      INTEGER i , j , k
      INTEGER klo , khi

      REAL plo , phi , tlo, thi , zlo , zhi
      REAL p_at_pconst , t_at_pconst , z_at_pconst
      REAL z_half_lowest

      LOGICAL  l1 , l2 , l3, found

      call get_wrf_p(cdfid,time)
      call get_wrf_z(cdfid,time)
      call get_wrf_qv(cdfid,time)
      call get_wrf_tk(cdfid,time)

!     Find least zeta level that is PCONST Pa above the surface.  
!     We later use this level to extrapolate a surface pressure 
!     and temperature, which is supposed to reduce the effect of 
!     the diurnal heating cycle in the pressure field.

      DO j = 1 , ny
         DO i = 1 , nx
            level(i,j) = -1
            k = 1
            found = .false.
            do while( (.not. found) .and. (k.le.nz))
               IF ( p(i,j,k) .LT. p(i,j,1)-PCONST ) THEN
                  level(i,j) = k
                  found = .true.
               END IF
               k = k+1
            END DO

            IF ( level(i,j) .EQ. -1 ) THEN
            PRINT '(A,I4,A)','Troubles finding level ', &
                  NINT(PCONST)/100,' above ground.'
            PRINT '(A,I4,A,I4,A)', &
                  'Problems first occur at (',i,',',j,')'
            PRINT '(A,F6.1,A)',    &
                  'Surface pressure = ',p(i,j,1)/100,' hPa.'
            STOP 'Error_in_finding_100_hPa_up'
         END IF

         END DO
      END DO

!     Get temperature PCONST Pa above surface.  Use this to extrapolate
!     the temperature at the surface and down to sea level.

      DO j = 1 , ny
         DO i = 1 , nx

            klo = MAX ( level(i,j) - 1 , 1      )
            khi = MIN ( klo + 1        , nz - 1 )

            IF ( klo .EQ. khi ) THEN
               PRINT '(A)','Trapping levels are weird.'
               PRINT '(A,I3,A,I3,A)','klo = ',klo,', khi = ',khi,   &
                     ': and they should not be equal.'
               STOP 'Error_trapping_levels'
            END IF

         plo = p(i,j,klo)
         phi = p(i,j,khi)
         tlo = tk(i,j,klo) + 0.608 * qv(i,j,klo)
         thi = tk(i,j,khi) + 0.608 * qv(i,j,khi)
!         zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j)
!         zhi = zetahalf(khi)/ztop*(ztop-terrain(i,j))+terrain(i,j)
         zlo = z(i,j,klo)
         zhi = z(i,j,khi)

         p_at_pconst = p(i,j,1) - pconst
         t_at_pconst = thi-(thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)
         z_at_pconst = zhi-(zhi-zlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)

         t_surf(i,j) = t_at_pconst*(p(i,j,1)/p_at_pconst)**(gamma*R/g)
         tsml(i,j) = t_at_pconst-gamma*z_at_pconst

         END DO
      END DO

!     If we follow a traditional computation, there is a correction 
!     to the sea level temperature if both the surface and sea level 
!     temnperatures are *too* hot.

      IF ( ridiculous_mm5_test ) THEN
         DO j = 1 , ny
            DO i = 1 , nx
               l1 = tsml(i,j) .LT. TC
               l2 = t_surf(i,j) .LE. TC
               l3 = .NOT. l1
               IF ( l2 .AND. l3 ) THEN
                  tsml(i,j) = TC
               ELSE
                  tsml(i,j) = TC - 0.005*(t_surf(i,j)-TC)**2
               END IF
            END DO
         END DO
      END IF

!     The grand finale: ta da!

      DO j = 1 , ny
        DO i = 1 , nx
!        z_half_lowest=zetahalf(1)/ztop*(ztop-terrain(i,j))+terrain(i,j)
          z_half_lowest=z(i,j,1)
          psml(i,j) = p(i,j,1) * EXP((2.*g*z_half_lowest)/        &
                                   (R*(tsml(i,j)+t_surf(i,j))))
        END DO
      END DO
      if (debug .GE. 500) then 
        print *,'sea pres input at weird location i=20,j=1,k=1'
        print *,'t=',tk(20,1,1),tk(20,2,1),tk(20,3,1)
        print *,'z=',z(20,1,1),z(20,2,1),z(20,3,1)
        print *,'p=',p(20,1,1),p(20,2,1),p(20,3,1)
        print *,'slp=',psml(20,1),psml(20,2),psml(20,3)
      endif 

      END SUBROUTINE compute_seaprs

      SUBROUTINE compute_rh(cdfid,time,var)

      IMPLICIT NONE
      INTEGER,intent(in)  :: cdfid, time 
      REAL,   intent(out) :: var(nx,ny,nz)

      REAL    svp1, svp2, svp3, svpt0
      parameter (SVP1=0.6112,SVP2=17.67,SVP3=29.65,SVPT0=273.15)
      integer i,j,k
      real qvs, es, pressure, temperature
      REAL    ep_2, r_d, r_v
      PARAMETER (r_d=287.,r_v=461.6, EP_2=R_d/R_v)
      REAL    ep_3
      PARAMETER (ep_3=0.622)
     
      call get_wrf_p(cdfid,time)
      call get_wrf_tk(cdfid,time)
      call get_wrf_qv(cdfid,time)

      do k=1,nz
        do j=1,ny
          do i=1,nx
!           es  = 1000.*svp1*
            es  = 10.*svp1*                         &
                  exp(svp2*(tk(i,j,k)-svpt0)/(tk(i,j,k)-svp3))
!           qvs = ep_2*es/(p(i,j,k)-es)
            qvs = ep_3*es/(0.01 * p(i,j,k)-(1.-ep_3)*es)
!        var = 100*amax1(1., qv(i,j,k)/qvs)
!       var(i,j,k) = 100.*qv(i,j,k)/qvs
            var(i,j,k) = 100.*AMAX1(AMIN1(qv(i,j,k)/qvs,1.0),0.0)
          enddo
        enddo
      enddo

      return
      END SUBROUTINE compute_rh

      SUBROUTINE compute_td2(cdfid,time,var)

      IMPLICIT NONE
      INTEGER,intent(in)  :: cdfid, time 
      REAL,   intent(out) :: var(nx,ny)

      REAL                :: q2(nx,ny),ps(nx,ny)

      ! compute q2, ps, then td
      if(debug .GE. 100) write(6,*)' get_wrf_theta: getting data for q2'
      call get_var_2d_real_cdf( cdfid, "Q2" , q2, nx, ny, time, debug)

      if(debug .GE. 100) write(6,*)' get_wrf_theta: getting data for Psfc'
      call get_var_2d_real_cdf( cdfid, "PSFC" , ps, nx, ny, time, debug)

      var = q2 * 0.01*ps / ( 0.622+q2) ! vapor pressure
      where (var < 0.01) var = 0.01    ! avoid problems near zero
 
      ! var = (243.5/( (17.67/log(var/6.112)) - 1.0)) ! Bolton's
                                                      ! approximation
      var = (243.5*log(var)-440.8)/(19.48-log(var))

      return

      END SUBROUTINE compute_td2

      SUBROUTINE compute_rh2(cdfid,time,var)

      IMPLICIT NONE
      INTEGER,intent(in)  :: cdfid, time 
      REAL,   intent(out) :: var(nx,ny)

      REAL                :: t2(nx,ny),q2(nx,ny),ps(nx,ny)

      REAL    svp1, svp2, svp3, svpt0
      parameter (SVP1=0.6112,SVP2=17.67,SVP3=29.65,SVPT0=273.15)
      integer i,j,k
      real qvs, es, pressure, temperature
      REAL    ep_2, r_d, r_v
      PARAMETER (r_d=287.,r_v=461.6, EP_2=R_d/R_v)
      REAL    ep_3
      PARAMETER (ep_3=0.622)
     
      if(debug .GE. 100) write(6,*)' get_wrf_theta: getting data for Psfc'
      call get_var_2d_real_cdf( cdfid, "PSFC" , ps, nx, ny, time, debug)

      if(debug .GE. 100) write(6,*)' get_wrf_theta: getting data for T2'
      call get_var_2d_real_cdf( cdfid, "T2" , t2, nx, ny, time, debug)

      if(debug .GE. 100) write(6,*)' get_wrf_theta: getting data for q2'
      call get_var_2d_real_cdf( cdfid, "Q2" , q2, nx, ny, time, debug)

        do j=1,ny
          do i=1,nx
            es  = 10.*svp1*                         &
                  exp(svp2*(t2(i,j)-svpt0)/(t2(i,j)-svp3))
            qvs = ep_3*es/(0.01 * ps(i,j)-(1.-ep_3)*es)
            var(i,j) = 100.*AMAX1(AMIN1(q2(i,j)/qvs,1.0),0.0)
          enddo
        enddo

      return
      END SUBROUTINE compute_rh2

END MODULE module_wrf_getvar
